% This file analyzes the experiments intended to show the phase transition
% [r,0,-r]

% here we measure the 2-form F, and combine with symmetry arguments to
% directly obtain the tensor gauge potential, and then the 3-form curvature

%% read all relevant data
% here we disable singular matrix warning during fitting
load('./data/FigS22.mat');
name_PT={'035','05','065','09','12','sqrt2','17'};
flag = {[3,1],[3,-1],[1,2],[1,-2]};
r=2e6;

for hlp_PT=1:7
    beta=0;phi=0;
    
    %% H_berry_curvature
    % calculate Berry curvature
    H_berry_curvature=-sum(F_munu_PT{hlp_PT},1)/2;
    H_berry_err=sqrt(F_munu_err_PT{hlp_PT}(1,:).^2+F_munu_err_PT{hlp_PT}(2,:).^2)/2;
    
    % Topological charge
    [QT,QT_err] = Herror_charge(alpha_sweep{hlp_PT},H_berry_curvature, H_berry_err);
    
    disp(['Q_T=',num2str(QT),' +- ',num2str(QT_err)])
    
    %% preparation for plot
    % calculate theore values
    A_pert=1e-7;
    alpha_fit=linspace(0,90,91)/180*pi;
    rabi_estimate_save = berry_rabi_freq(alpha_fit,beta,phi,h(hlp_PT),r,flag);
    rabi_estimate_save=rabi_estimate_save/1e6;
    
    F_munu_theory=zeros(3,length(alpha_fit));
    
    for hlp=1:length(alpha_fit)
        alpha=alpha_fit(hlp);
        F_munu_theory(:,hlp)=F_munu_solve(alpha,beta,phi,h(hlp_PT),A_pert);
    end
    
    F_munu_theory=real(F_munu_theory(2:3,:));
    
    H_munu_theory=-real(sum(F_munu_theory,1)/2);
    QT_theory = Herror_charge(alpha_fit,H_munu_theory, zeros(size(H_munu_theory)));
    
    %% plot
    % Rabi oscillations
    figure
    cc=colormap(piratepal(9,1));
    cc=cc([1,2,3,6],:);
    
    % figure
    % SQ
    hold on
    for hlp=1:length(flag)
        plot(alpha_fit/pi*180,rabi_estimate_save(hlp,:,1),'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
    end
    for hlp=1:length(flag)
        errorbar(alpha_sweep{hlp_PT}/pi*180,rabi_fit_save_PT{hlp_PT}(hlp,:,1)./E_pert_save_PT{hlp_PT}(hlp,:,1),drabi_fit_save_PT{hlp_PT}(hlp,:,1)./E_pert_save_PT{hlp_PT}(hlp,:,1),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
    end
    
    ymax=(max(max(rabi_estimate_save(:,:,1))))+0.25;
    ylim([0,ymax])
    xlim([0,90])
    xticks(0:15:90)
    xlabel('\alpha(\circ)','FontSize',18)
    ylabel('Matrix Element (MHz)','FontSize',18)
    title('SQ','FontSize',20)
    hold off
    set(gca,'FontSize',18)
    
    lgd=legend({'$$\Gamma_{\phi\alpha}$$','$$\Gamma_{\phi\bar{\alpha}}$$','$$\Gamma_{\alpha\beta}$$','$$\Gamma_{\alpha\bar{\beta}}$$'});
    set(lgd,'Interpreter','latex');
    lgd.FontSize=18;
    set(lgd,'Location','best')
    
%     saveas(gcf,['F_SQ_',name_PT{hlp_PT}],'epsc')
    % DQ
    figure
    hold on
    for hlp=1:length(flag)
        plot(alpha_fit/pi*180,rabi_estimate_save(hlp,:,2),'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
    end
    for hlp=1:length(flag)
        errorbar(alpha_sweep{hlp_PT}/pi*180,rabi_fit_save_PT{hlp_PT}(hlp,:,2)./E_pert_save_PT{hlp_PT}(hlp,:,2),drabi_fit_save_PT{hlp_PT}(hlp,:,2)./E_pert_save_PT{hlp_PT}(hlp,:,2),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
    end
    
    ymax=(max(max(rabi_estimate_save(:,:,2))))+0.25;
    ylim([0,ymax])
    xlim([0,90])
    xticks(0:15:90)
    xlabel('\alpha(\circ)','FontSize',18)
    title('DQ','FontSize',20)
    hold off
    set(gca,'FontSize',18)
    
    title('DQ','FontSize',20)
    hold off
    set(gca,'FontSize',18)
%     saveas(gcf,['F_DQ_',name_PT{hlp_PT}],'epsc')
    
    % 2-form and 3-form Berry curvature
    figure
    cc_g=colormap(piratepal(9,1));
    hold on
    for hlp=1:2
        plot(alpha_fit/pi*180,F_munu_theory(hlp,:),'LineWidth',1.5,'Color',cc_g(hlp,:),'HandleVisibility','off')
    end
    for hlp=1:2
        errorbar(alpha_sweep{hlp_PT}/pi*180,F_munu_PT{hlp_PT}(hlp,:),F_munu_err_PT{hlp_PT}(hlp,:),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc_g(hlp,:))
    end
    hold off
    xlabel('\alpha(\circ)','FontSize',18)
    title('F_{\mu\nu}','FontSize',20)
    set(gca,'FontSize',18)
    xlim([0,90])
    xticks(0:15:90)
    % ylim([0,6])
    
    lgd=legend({'$\mathcal{F}_{\phi\alpha}$','$\mathcal{F}_{\alpha\beta}$'},'Interpreter','latex');
    lgd.FontSize=18;
    
%     saveas(gcf,['F_',name_PT{hlp_PT}],'epsc')
end


function [rabi_g,E_total] = berry_rabi_freq(alpha_sweep,beta,phi,h,r,flag)

rabi_g = zeros(length(flag),length(alpha_sweep),2);
E_total = zeros(2,length(alpha_sweep)); % 1: SQ, 2: DQ

for hlp=1:length(alpha_sweep)
    alpha=alpha_sweep(hlp);
    [V,E]=Ham_detune(alpha,beta,phi,h);
    um=V(:,1);u0=V(:,2);up=V(:,3);
    
    E_total(:,hlp)=[abs(E(2)-E(1))*2;abs(E(3)-E(1))]*r; % note we already take 2X the SQ frequency here
    
    for hlp_flg=1:length(flag)
        flg_tmp=flag{hlp_flg};
        
        for hlp_omega=1:2 % omega
            dH1 = DHam(alpha,beta,phi,flg_tmp(1));
            dH2 = DHam(alpha,beta,phi,flg_tmp(end));
            
            dH=dH1+1i*dH2;
            
            switch hlp_omega
                case 1
                    rabi_g(hlp_flg,hlp,hlp_omega)=abs(um'*dH*u0*r);
                    
                case 2
                    rabi_g(hlp_flg,hlp,hlp_omega)=abs(um'*dH*up*r);
            end
        end
    end
end
end

function Ham_out = DHam(alpha,beta,phi,D_order) % calculates the derivative
switch abs(D_order)
    case 1
        Ham_out=[0,-sin(alpha)*exp(-1i*beta),0;-sin(alpha)*exp(1i*beta),0,cos(alpha)*exp(1i*phi);...
            0,cos(alpha)*exp(-1i*phi),0];
    case 2
        Ham_out=[0,-1i*cos(alpha)*exp(-1i*beta),0;1i*cos(alpha)*exp(1i*beta),0,0;...
            0,0,0];
    case 3
        Ham_out=[0,0,0;0,0,1i*sin(alpha)*exp(1i*phi);...
            0,-1i*sin(alpha)*exp(-1i*phi),0];
end

if D_order<0
    Ham_out=-Ham_out; % mu,-nu
end
end

function B_munu = B_munu_solve(alpha,beta,phi,h,A_pert)
% numerically solve for B=phi F
[F_munu,vg] = F_munu_solve(alpha,beta,phi,h,A_pert);

B_munu=F_munu*(-1i/2)*log(vg(1)*vg(2)*vg(3));
end

function [F_munu,vg] = F_munu_solve(alpha,beta,phi,h,A_pert)
% numerically solve for Berry curvature F_munu
F_munu=zeros(3,1);
flag = {[2,3],[3,1],[1,2]};

[A_mu,vg] = A_mu_solve(alpha,beta,phi,h,A_pert);

for hlp_flg=1:length(flag)
    flg_tmp=flag{hlp_flg};
    
    % \partial_mu A_\nu
    param_mod=[0,0,0];param_mod(flg_tmp(1))=A_pert;
    alpha_pert=alpha+param_mod(1);beta_pert=beta+param_mod(2);phi_pert=phi+param_mod(3);
    
    A_mu_pert = A_mu_solve(alpha_pert,beta_pert,phi_pert,h,A_pert);
    
    F_munu(hlp_flg)=(A_mu_pert(flg_tmp(2))-A_mu(flg_tmp(2)))/A_pert;
    
    % - \partial_nu A_\mu
    param_mod=[0,0,0];param_mod(flg_tmp(2))=A_pert;
    alpha_pert=alpha+param_mod(1);beta_pert=beta+param_mod(2);phi_pert=phi+param_mod(3);
    
    A_mu_pert = A_mu_solve(alpha_pert,beta_pert,phi_pert,h,A_pert);
    
    F_munu(hlp_flg)=F_munu(hlp_flg)-(A_mu_pert(flg_tmp(1))-A_mu(flg_tmp(1)))/A_pert;
end
end

function [A_mu,vg] = A_mu_solve(alpha,beta,phi,h,A_pert)
% numerically solve for Berry connection A_mu
A_mu = zeros(3,1);
[uv,~]=Ham_detune(alpha,beta,phi,h); % eigenvectors
vg=uv(:,1);

for k=1:3
    param_mod=[0,0,0];param_mod(k)=A_pert;
    
    [uv_pert,~] = Ham_detune(alpha,beta,phi,h,param_mod);
    vg_pert = uv_pert(:,1);
    
    A_mu(k)=1i*vg'*(vg_pert-vg)/A_pert;
end
end

function [Vout,Eout] = Ham_detune(alpha,beta,phi,h,A_pert) %
if nargin>4
    alpha=alpha+A_pert(1);beta=beta+A_pert(2);phi=phi+A_pert(3);
end

Ham_out=[0,cos(alpha)*exp(-1i*beta),0;cos(alpha)*exp(1i*beta),0,sin(alpha)*exp(1i*phi);...
    0,sin(alpha)*exp(-1i*phi),0];
Ham_out=Ham_out+diag([h,0,-h]);


[V,E]=eig(Ham_out);
[Eout,idx]=sort(diag(E));
V=V(:,idx);
um=V(:,1)/sqrt(V(:,1)'*V(:,1));
u0=V(:,2)/sqrt(V(:,2)'*V(:,2));
up=V(:,3)/sqrt(V(:,3)'*V(:,3));

if abs(um(2))>1e-5
    um=um/um(2)*abs(um(2));
end
if abs(u0(2))>1e-5
    u0=u0/u0(2)*abs(u0(2));
end
if abs(up(2))>1e-5
    up=up/up(2)*abs(up(2));
end

Vout=[um,u0,up];

end



